function [hs,e,j] = calcHs_TDIP(invP,etaModel,useMumps)

k = 1/invP.disT(1,2);

M = k*invP.mu*speye(invP.mesh.ne);
S = sdiag(invP.mesh.Af*(invP.sigmaModel.*(1 + k*invP.tauModel)));
St = sdiag(invP.mesh.Af*(k*invP.tauModel.*invP.sigmaModel));
T = sdiag(invP.mesh.Af*(1 + k*invP.tauModel.*(1 - etaModel)));
Tt = sdiag(invP.mesh.Af*(k*invP.tauModel.*(1-etaModel)));

A = invP.mesh.CURL'*(S\T)*invP.mesh.CURL + M;


hs = cell(invP.nt,1);
e = cell(invP.nt,1);
j = cell(invP.nt,1);

prev_hs = invP.h0;
prev_e = zeros(invP.mesh.nf,invP.ntx);
prev_j = zeros(invP.mesh.nf,invP.ntx);

for ii = 1:invP.nt
    fprintf(1,'Solving time %i of %i....\t',ii,invP.nt);
    
    rhs = k*invP.mu*prev_hs - invP.mesh.CURL'*(S\St)*prev_e ...
        + invP.mesh.CURL'*(S\Tt)*prev_j;
    
    hs{ii} = A\rhs;
    e{ii} = (S\T)*invP.mesh.CURL*hs{ii} + (S\St)*prev_e - (S\Tt)*prev_j;
    j{ii} = (T\S)*e{ii} - (T\St)*prev_e + (T\Tt)*prev_j;
    
    prev_hs = hs{ii};
    prev_e = e{ii};
    prev_j = j{ii};
    
    fprintf(1,'Done\n');
end
